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The non-linear transformations incurred by the rays in an optical system can be suit- 
ably described by matrices to any desired order of approximation. In systems composed 
of uniform refractive index elements, each individual ray refraction or translation has an 
associated matrix and a succession of transformations correspond to the product of the re- 
spective matrices. This paper describes a general method to find the matrix coefficients for 
translation and surface refraction irrespective of the surface shape or the order of approxi- 
mation. The choice of coordinates is unusual as the orientation of the ray is characterised 
by the direction cosines, rather than slopes; this is shown to greatly simplify and generalise 
coefficient calculation. Two examples are shown in order to demonstrate the power of the 
method: The first is the determination of seventh order coefficients for spherical surfaces 
and the second is the determination of third order coefficients for a toroidal surface. 
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1 Introduction 

An optical system can be effectively modelled in paraxial approximation by a product of 4 x 4 matrices, 
each representing one elementary transformation of the light rays [ 1 1 ; the elementary transformations 
are either translations of the ray in homogeneous media or the effects of surfaces separating different 
media. A more accurate approach implies the consideration of higher order terms but the fact that 
Snell's law makes use of the sine function rules out the terms of even order; as a consequence, when 
one wants to improve on the paraxial approximation, one has to consider third order terms. 

Aberrations have already been studied extensively J2][3) but work is still going on in order to design 
symbolic models of optical systems that computers can use for optimisation purposes and humans can 
look at to gain a better understanding of systems' performance. 
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The matrix theory has been extended to deal with higher order terms [4| through the use of a vector 
basis that incorporates two position and two orientation coordinates as well as all their third or higher 
order monomials, increasing the overall dimension which becomes 24 for third order approximation. 
It is possible to apply axis symmetry to reduce the matrix dimension through the use of complex 
coordinates and their higher order monomials; for instance, in third order the matrices to be considered 
for axis-symmetric systems are 8x8 0|5]|. 

The set of four coordinates normally used to describe the ray consists of the two rectangular coordi- 
nates, x and y along with the two ray slopes u = dx/dz and v = dy/dz. I chose to replace the ray 
slopes with the direction cosines relative to the coordinate axes, respectively s and t, in order to allow 
an easier and more elegant formulation of the Snell's law at a surface but at the expense of rendering 
the translation transformation non-linear. 

This work details a general method that can be used to establish all the coefficients needed for any 
order modelling of optical systems built with surfaces of unspecified shape. The power of the method 
is exemplified with determination of matrix coefficients for spherical systems in seventh order and for 
a toroidal surface in third order. 



2 The choice of coordinates 

The optical axis is assumed to lie along the z axis, so the position of any point is determined by the two 
coordinates x and y; when dealing with axis symmetric systems, the two coordinates can be combined 
in the complex coordinate X = x + iy. The ray orientation is defined by cosines of the angles with the 
x and y axes, respectively s and t; again a complex orientation coordinate can be used, S = s + it, in 
order to simplify the rotational treatment of axis symmetric systems. 

Using the same notation of reference |4|, x& is the vector whose components are the coordinates 
and their monomials up to order n: 
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where total order j ' + k + / + m is n or less and odd, with all exponents greater than zero. By convention 
the vector elements are placed in the order of smaller n and larger four digit number jklm. 

Any transformation of x& into x'& can be represented by a square matrix S of dimension N equal 
to the size of x&, following the equation: 



x'& = Sx& 



(2) 
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The matrix S has the following structure: 

a — ( f*4x4 H 4x ( N _4) 

V 0(N-4)x4 E( N _4) X ( N _4) 

The four submatrices that form S have the following meaning: P4 X 4 contains the paraxial constants 
of the optical system, H 4x ( N _ 4 ) has A(N — 4) high order coefficients from the series expansion of the 
transformed coordinates, 0( N _ 4 ) x4 is composed of zeroes and finally E(n-4)x(N-4) is obtained from 
the elements of the first four lines by a procedure called extension whereby the elements of lines 5 to 
N from x'& are calculated and all terms of orders higher than n are dropped. More specifically, if we 
wish to determine the coefficients for the line corresponding to the monomial x^y k s l t m , we take the 
polynomial expressions for x, y, s and t from lines 1 to 4 and raise them to the powers j, k, I and m, 
respectively, making their product afterwards; the result is a polynomial of degree n x (j : + k + I + m) 
from which only the terms of orders up to n should be considered. The submatrix E( N _ 4 ) X ( N _ 4 ) can 
itself be subdivided into components of different orders and components with all zero elements. 

Although the ray transformation is described by an ./V x N matrix, the above considerations show 
that only the A(N — 4) elements from the first four lines need be considered, as the extension procedure 
is standard for all transformations. In cases where symmetry exists the number of independent coeffi- 
cients can be greatly reduced 0|6). It is apparent that the matrix for any given transformation can be 
considered completely defined when the first four lines' coefficients have been evaluated, i.e. when the 
transformed coordinates have been expressed in a power series of order n. 

The coordinate conventions made above are not the same as those made in references [4 5 1 and 
indeed by the majority of authors; they may appear to be a poorer choice because, as we shall see, they 
lead to a non-linear translation transformation; on the other hand they will simplify the determination 
of refraction coefficients. It should be noted, however, that a coordinate change can be seen as a 
transformation, itself governed by a matrix; it is not difficult to follow the procedure outlined below 
and then change from cosines into slopes, if needed. 




3 The elementary transformations 
3.1 The translation matrix 

The first elementary transformation that has to be considered is just the translation of the ray in an 
homogeneous medium; the orientation coordinates don't change but the position coordinates change 
according to the equations 

, se 
x = x + 



Vl - s 2 - t 2 ' 
vl "«> < 2 (4) 



where e is the distance travelled, measured along the optical axis. 
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The series expansion of the equations is rather straightforward; up to the seventh order it is: 

j €■ q (3 n 3(3 c 6(3 q o 3c /I 

x' = x + es + -s d + -st 2 + — s b + — s'H 2 + — sF 
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The previous equations give directly the coefficients for lines 1 and 2 of the matrix T which describes 
the translation transformation; lines 3 and 4 are made up of zeroes, except for the diagonal elements 
which are 1; this translates the two relations s' = s and t' = t. The lines 5 to N can be obtained by the 
extension procedure described before. 

If slopes were used instead of direction cosines the translation matrix would have a much simpler 
form, corresponding to a linear transformation [4]. 

3.2 Refraction formulae 

The influence of the surface power on the ray coordinates must now be considered. It is useful to 
analyse separately the changes on the ray orientation and the modification of the position coordinates. 
The ray orientation changes due to the application of Snell's law when the ray encounters the surface 
of separation between the two media; let this surface be defined by the general formula: 

i(x,y,z) = (6) 

Let a ray impinge on the surface on a point of coordinates x and y with direction cosines s and t. 
The normal to the surface at the incidence point is the vector: 

n = grad f ; (7) 

the direction of the incident ray is the unit vector: 

s 

t I ; (8) 
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and the refracted ray has a direction which is represented by the vector: 



v' 



t' , (9) 



VI - s' 2 - t' 2 



where s' and t' are the direction cosines after refraction. 

One can apply Snell's law by equating the cross products of the ray directions with the normal on 
both sides of the surface multiplied by the respective refractive indices: 

i^v§<)n = v / (g | n, (10) 

where u represents the refractive index ratio of the two media. 
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The algebraic solutions of the equation above can be found for many surfaces by means of suitable 
symbolic processing software like Mathematica; in certain cases the solution is easier to find after a 
suitable change of coordinates along the z axis, which will not affect the final result; for instance, for 
a spherical surface it is convenient to shift the coordinate origin to the center. In order to find the 
coefficients for the matrix Ri , representing the change in the ray orientation, the explicit expressions 
for s' and t' must be expanded in series; this can again be programmed into Mathematica. The examples 
given later show this procedure for two different surfaces. 



3.3 Surface offset 

The use of just two position coordinates implies that these are referenced to a plane normal to the 
axis. When a ray is refracted at a surface, the position coordinates that apply are those of the incidence 
point; nevertheless the translation to the surface is referenced to the vertex plane and thus the position 
coordinates that come out of the translation transformation, are those of the intersection with that plane; 
the difference between the two positions is named surface offset. 



Vertex plane 



Surface 




Optical axis 



Figure 1 : The ray intersects the surface at a point X\ which is different both from the point of inter- 
section of the incident ray with the plane of the vertex, X, and the point of intersection of 
the refracted ray with the same plane, X?,. The surface is responsible for three successive 
transformations: 1 - an offset from X to X\, 2 - the refraction and 3 - the offset from X\ to 
X 2 . 

The figure illustrates the problem that must be solved: The reference plane for the position coor- 
dinates of the ray in a refraction process is the plane of the surface vertex; the coordinates of the point 
of intersection of the ray with this plane are not the same before and after the refraction. The surface 
matrix must account not only for the orientation changes but also for the position offset introduced by 
the refraction process. 

If the origin of the rectangular coordinates is located at the vertex and (x, y, 0) are the coordinates 
of the ray when it crosses the plane of the vertex, the current coordinates of a point on the incident ray 
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are given by the two equations: 



x' = x + 



VI - s 2 - t 2 ' 



As the intersection point has to verify both dTH and © it is possible to solve a system of three 
simultaneous equations to determine its coordinates. In fact there may be several intersection points, 
as the ray may intercept a complex surface in several points; it is not difficult, however, to select 
the solution of interest by careful examination of the surface region that has been hit; we are usually 
interested in the solution that has the smallest \z'\ value. 

As before a series expansion of the exact expressions of x' and y' must be performed taking only 
the terms up to the nth order; this will yield the coefficients for the matrix Ti which describes the 
transformation from point X to point X\ and will be called the forward offset. 

The offset from point X\ to point X2, designated backward offset, is the reverse transformation of 
the forward offset, for a ray whose direction coordinates have been modified by refraction. We solve 
(fTTT> and ® i n terms of x and y and apply the previous procedure to determine the coefficients of the 
transformation matrix, T2. 

The matrix describing the transformations imposed by the surface is the product T2R1T1, as the 
ray has to suffer the three successive transformations: forward offset, refraction and backward offset. 
The surface matrix is called R. 



4 Stop effects 

The entrance pupil of the optical system plays an important role in the overall aberrations; the center of 
the pupil defines the chief ray, which determines the paraxial image point. The ray fan usually extends 
equally in all directions around the chief ray and the points of intersection of the various rays in the 
fan with the image plane determine the ray aberrations. When matrices are used to model the system 
the image appears described in terms of the position and orientation coordinates of the rays when these 
are subjected to the first transformation, be it a refraction or a translation; in terms of ray aberrations 
this corresponds to a situation where the object is located at infinity and the entrance pupil is placed 
where the first transformation takes place. In fact, the orientation coordinates play the role of object 
coordinates and the position coordinates are the actual pupil coordinates, if this coincides with the first 
transformation; the chief ray can be found easily just by zeroing the position coordinate. 

Ray aberrations can be evaluated correctly if the image is described in terms of both object and stop 
coordinates or some appropriate substitutes; an object point is then set by the object coordinates, the 
paraxial image point is found when the stop coordinate is zero and the aberrations are described by the 
differences to this point. This can be done adequately when the stop is located before the first surface 
and so coincides with the entrance pupil. If the object is at infinity the problem is very easy to solve and 
application examples have already been described [6|. If the entrance pupil is located at the position z s 
relative to the first surface, this means that the rays must perform a translation of — z s prior to hitting it; 
this can be accommodated via the right product by a translation matrix T_ Zs , where the index indicates 
the amount of translation. 

The case of near objects is more difficult to deal with because we are faced with a dilemma: If a 
translation from the pupil position to the first surface is applied the position coordinates of the incident 
rays become pupil coordinates but the rays' origins on the object are lost. Conversely we can apply a 
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Figure 2: The coordinates of the point where the ray crosses the entrance pupil, together with the object 
point coordinates are a convenient way of defining the ray orientation in order to incorporate 
stop effects. 



translation from the object position, ending up with the reverse problem of knowing the rays' origins 
and ignoring their point of passage through the pupil. We need the position coordinates both at the 
object and pupil location; this suggests that the ray's orientation should be specified not by its direction 
cosines but by the coordinates of the point of passage through the pupil. In figure|2]a ray leaves an object 
point of coordinates (x, y) and crosses the entrance pupil plane at a point of coordinates (x p , y p ); then 
the direction cosines, s and t can be calculated by: 



rp rp *\ 2 I *y 2 



y P -y 



\J {vp - y) 2 + 4 



(12) 



where z p is the position of the pupil relative to the object. 

The matrix theory developed before was based on direction cosines rather than pupil coordinates; in 
consequence a coordinate change is needed before the translation from the object to the first surface is 
applied to the rays. Rather than use equations (fT2l it is more convenient to look at them as a coordi- 
nate change governed by a matrix T p ; for the determination of its coefficients we apply the standard 
procedure of series expansion. The following is the result for 7th order: 

1 1 1 o 3 o . 3 o 3 ^ 15 4 

q rp I rp I *-J rp rp I rp rp ™ rp rp I rp rp 
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15 q o 15 o q 15 a 3 c 5 7 35 c 105 c 9 

rp^rp^ _|_ rp ^ rp ° rp rp^ _|_ rp ° _|_ rp ' ™ U rp _|_ rp° rp^ 

4zJ p + 4z| X X p 8z| XX p + 8z| p + 164 16 4 16 4 P 

175 4 o 175 o 4 105 o * 35 fi 5 7 

/y - y ■ T I v ' ' rp rp rp I rp rp rp ' I 1 T I 

16zf p + 164 p Wz f p 16 4 p 1Q 4 p ' 

a similar expression exists for t. The equation above allows the determination of the coefficients for 
matrix T p which converts the coordinates (x, y, x p , y p ) into (x, y, s, t) before the rays enter the system 
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at the object position; as a symmetrical conversion is not performed when the rays leave the system, 
there is no place for the use of matrix Tp 1 , as would be the case in a similarity transformation. 
A case study would involve the following successive steps: 

• Determine the matrices for all the translation and refraction transformations suffered by the ray 
from the object point to the image plane. 

• Determine the matrix T p , relative to the entrance pupil position. 

• Multiply all the matrices in reverse order, starting with the translation from the last surface to the 
image plane and ending with matrix T p ; the result is matrix S. 

• Determine a vector base x& with the set of coordinates (x, y, x p , y p ). 

• Make the product x'& = Sx&. 

Once x'& has been calculated, its first four elements are polynomials of nth degree on the indepen- 
dent variables {x, y, x p , y p ), which model the position and direction cosines of the rays on the image 
plane. 

Most optical systems have stops located after the first surface, and so the method described above 
is not adequate. One can always find the gaussian entrance pupil of any system and so revert to the 
previous situation; this will work even if the entrance pupil is found to lie beyond the first surface, in 
which case it will give rise to a negative translation. The problem is that the gaussian entrance pupil 
is a first order transformation of the stop, whereas the real entrance pupil is an aberrated image of the 
stop given by the portion of the optical system that lies before it. Ideally one should try to describe the 
image in terms of the object coordinates and the point of passage through the stop, which can be rather 
difficult. Alternatively one can find a reference ray, defined as the ray that passes through the center 
of the stop and is not necessarily the same as the chief ray, and work with position and orientation 
coordinates on the image plane, relative to the reference ray. The plot of the relative ray position versus 
its relative direction cosines can be used to describe ray aberrations Q. 

5 Coefficients for a spherical surface 

A sphere of radius r with its center at z = r is defined by the equation: 

x 2 + y 2 + z 2 - 2zr = 0; (14) 

this equation will replace © in the algorithm for the determination of matrix coefficients. In a previous 
work [ 6 1 we presented the results for third order; these are now extended to seventh order. 

This is a situation where symmetry must be considered in order to reduce the size of the matrices 
involved; we use the complex variables X and S previously defined, leading to a complex vector base 
X& whose elements have the general form X J X* k S l S* m with X* and S* representing the complex 
conjugates of X and S. Kondo et al. [4| have shown that the powers must obey the condition j — k + 
I — m = 1; for seventh order the allowable combinations are: 

1000, 0010, 2100, 2001, 1110, 1011, 0120, 0021, 3200, 3101, 3002, 2210, 2111, 2012, 1220, 1121, 
1022, 0230, 0131, 0032, 4300, 4201, 4102, 4003, 3310, 3211, 3112, 3013, 2320, 2221, 2122, 2023, 
1330, 1231, 1132, 1033, 0340, 0241, 0142, 0043 . 
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Equation (flOl can now be solved and the resulting expressions for s' and t' can be expanded in series 
up to the seventh order; the results can be combined into one complex variable S'. The second column 
of table^hsts the coefficients for this expansion. Similarly the forward offset coefficients can be found 
solving equations dTH and combining x' and y' expansions into one single complex variable X\ ; the 
resulting coefficients are listed in column 3 of tabled The backward offset coefficients are obtained in 
a similar way, although now it is the x and y expansions that must be combined into a single complex 
variable X2. 



Toroidal surfaces are used to correct eye astigmatism; they can be fabricated easily without resorting 
to sophisticated machinery. In this example we deem to show that the method is applicable to general 
surfaces but it is not our purpose to list high order coefficients for this particular shape; therefore we 
will restrict the study to third order. Also, as there is no rotational symmetry we will not use complex 
coordinates; there is symmetry relative to two perpendicular planes and if we were to use complex 
coordinates the terms should obey the condition j — k + l — to = ±1 as explained in reference 0; we 
feel that it is clearer not to invoke symmetry at all. 

The surface is generated by a circle of radius r\ normal to the xz plane whose center is displaced 
along a circle of radius T2 laying on this plane. This surface has curvature radii of r\ and r\ + r2 
respectively on the xz and yz planes. The surface equation can be written as: 



The solution of equation (flOl originates the expressions for s' and t'; these are then expanded in 
series up to the third order and the resulting coefficients are listed in columns 2 and 3 of table |2] The 
same procedure was applied to the expressions for x' and y' resulting from equation (fTTT i: columns 4 
and 5 of the same table list these coefficients. The backward offset coefficients need not be listed, as in 
third order they can be obtained by sign reversal from the forward offset coefficients. 

7 Discussion and conclusion 

Matrices can be used to model optical systems built with surfaces of unspecified shape to any desired 
degree of approximation. The method that was described differs from those found in literature by the 
choice of direction cosines rather than ray slopes to specify ray orientations. This results in a very ele- 
gant formulation of Snell's law which can be easily programmed into symbolic computation software. 
In spite of some added complication in translation matrices the algorithm is limited only by computing 
power in it's ability to determine the matrix coefficients, irrespective of surface complexity and order 
of approximation. For axis-symmetric systems the matrices have a dimension 40 x 40 for seventh 
order and can usually be dealt with by ordinary PCs; the influence of stops can also be incorporated 
in the matrix description as was demonstrated. The unconventional coordinate system that was used 
can, if needed, be converted to the more usual ray slope coordinate system through the product with 
conversion matrices of the same dimension. 

Two examples were shown for illustration of the method's capabilities; the first one supplies a list 
of seventh order coefficients for systems made up with spherical surfaces and allows the reader an 
immediate use. The second example applies the method to a fourth degree toroidal surface, showing 
that complex shapes can be accommodated. 



6 Coefficients for an asymmetric surface 




(15) 
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Table 1 : 7th order coefficients for spherical surfaces 



Monomial 


Refraction (S f ) 
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Backward offset (X2) 
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Table 2: 3rd order coefficients for a torous 



Monomial Refraction (V) Refraction (t') Forward offset (x') Forward offset (y 1 ) 
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